#ifndef _QGCM_SSYSTEM_
#define _QGCM_SSYSTEM_

#include <list>
#include <netcdfcpp.h>
#include "datatypes.h"
#include "constants.h"
#include "forward_classes.h"
#include "static.h"
#include <mpParser.h>
using namespace mup;
// -----------------------------------------------------------------------------------------
struct lopinfo {
  public:
    // - - - - - - - - - - - -
    struct element { 
      INDEX angular_mode;
      rank2 field_indices;
      ParserX * parser;
      // - - - - - - - - - - - -
      element();
      element(const element & ol);
      element & operator=(const element & ol);
      ~element();
      void SetExpr(const token_list & expr);
      // - - - - - - - - - - - -
      REAL imag() const;
      REAL real() const;
      void assign(const INDEX & angm, const rank2 & fpos, const token_list & expr, const ParserX & Pinit);
    }; // end element class 
    // - - - - - - - - - - - -
    vector< element > lop_list;
    lopinfo();
    ~lopinfo();
}; // end lopinfo class
// -----------------------------------------------------------------------------------------
struct linearOPtranslationdata {
	string OPname;
	INDEX a_mode;
	bool evaluate_op;    // are we to evaluate this operator
	vector< pair< rank3, token_list > > expressions; // rank3 is the position (m,i,j) 
};
// -----------------------------------------------------------------------------------------
struct initialdata {
	// - - - - - - - - - - - - - - - - - - - - - - - - 
	// variables beneath are iterative for populating
	// the fields
	// - - - - - - - - - - - - - - - - - - - - - - - -
	Value index_i;
	Value index_j;
	Value index_k;
	Value index_mode;
	Value modeval; // actual value of the mode
	Value kval; // floating point value K 

	struct rule {
		//enum dependency_type { field, coord, angmode };
		ParserX parser;
		//vector<dependency_type> order;
		rank4 rpos;
		// - - - - - - - - - - - - - - - - - - - - - - -  
		string namelist; // name of namelist where rule is defined
		int line_number; // where expression begun
		rule();
		~rule();
	};
	vector<rule> rules;
	initialdata();
	~initialdata();
	private:
	initialdata(const initialdata &); // no implicit copies allowed
	initialdata operator=(const initialdata &); // no implicit assignments allowed
};
// *****************************************************************************************
class spreprocdata {
	friend class ssystem;
	friend class model_class;
	friend class compiler;
	friend class debugger;
	// class stores preprocessing data needed for the compilation of the ssystems's
	// linear operators into an executable state 
	public:
	spreprocdata();
	~spreprocdata();
	initialdata * getdata();
	const vector<linearOPtranslationdata> & getOPdata() const;
	void clear();
	private:
	vector<linearOPtranslationdata> OPdata;
	vector< pair<INDEX, token_list> > dissip_data;
	struct ihomo_data {
		string name;
		rank3 pos;
		string ang_mode_var_name;
		string k_var_name;
		string field_i_var_name;
		string field_j_var_name;
		token_list expr;
	};
	vector< ihomo_data  > inhomogeneous_data;  
	initialdata * initial_conditions;
};

// *****************************************************************************************
// Class ssystem ...
// *****************************************************************************************
class ssystem {
	friend class psystem;
	friend class model_class;
	friend class compiler;
	public:

	void setstepequaltimeonly(bool); // function to set var stepequaltimeonly
	bool save_model_state(NcFile *) const;
	bool restore_model_state(NcFile *);

	private:
	class clocal;
	clocal * local;

	explicit ssystem(const ssystem &); // no copy constructor should be callled!
	ssystem & operator=(const ssystem &); // no equal assignments allowed

	public:
	REAL getelapsedtime(const unsigned int &) const;
	void setnumbertimesteps(const unsigned int &);
	void settimestepwidth(const unsigned int &, const REAL &);

	const REAL & gettimestep(const unsigned int &) const;
	void rewindtime(const unsigned int &);

	void stepforward(REAL timestep, bool initialguess);
	static_container * getRcoeffs_dyn();  // for modifying (NOT THREAD SAFE!!)
	static_container * getScoeffs_dyn();  // for modifying (NOT THREAD SAFE!!)

	void settprimeintegrationlist(const vector<int> *);

	REAL getflux(unsigned int index_i, unsigned int index_j, unsigned int tindex, bool initial_guess) const;

	public:
	ssystem();
	~ssystem();

	const qgcomplex & get_register(const regset&, const INDEX, const INDEX, const INDEX, const INDEX) const;

	const regset & get_R_set(unsigned int t, unsigned int t_primed) const;
	const regset & get_S_set(unsigned int t, unsigned int t_primed) const;
	const regset & get_phi_set(unsigned int t, unsigned int t_primed) const;
	const regset & get_G_set(unsigned int t, unsigned int t_primed) const;
	const regset & get_F_phi_set(unsigned int t, unsigned int t_prime) const;
	const regset & get_F_G_set(unsigned int t, unsigned int t_prime) const;
	const regset & get_F_phi_regularized_set(unsigned int t) const;
	const regset & get_R_regularized_set(unsigned int t) const;
	const regset & get_S_regularized_set(unsigned int t) const;

	const twotimevar & get_R_data() const;
	const twotimevar & get_S_data() const;
	const twotimevar & get_phi_data() const;
	const twotimevar & get_G_data() const;
	const twotimevar & get_F_phi_data() const;
	const twotimevar & get_F_G_data() const;
	const twotimevar & get_F_phi_regularized_data() const;
	const twotimevar & get_R_regularized_data() const;
	const twotimevar & get_S_regularized_data() const;

	const vector< pair<INDEX, ParserX> > & get_dissipation_ops() const;
        const vector< pair<rank3, ParserX> > & get_inhomogeneous_ops() const;
        const lopinfo & get_linear_op_info() const;

	bool is_model_using_regularization() const;
	int get_tindex() const;
	void sync_tindex(const int); // synchronize tindex with the P-system in case of INT signal

	void init(model_class & mo, spreprocdata &);
	const bool isinitialized() const;
	void fregdump(FILE *, const regset &) const;
	void dumpdata() const;

#ifdef DEBUG
	void regdump(const regset &) const;
	void dumpR(unsigned int t, unsigned int t2) const;
	void dumpS(unsigned int t, unsigned int t2) const;
	void dumpPhi(unsigned int t, unsigned int t2) const;
	void dumpG(unsigned int t, unsigned int t2) const;
	void dumpFPhi(unsigned int t, unsigned int t2) const;
	void dumpFG(unsigned int t, unsigned int t2) const;

	void dumpFPhi_regularized(unsigned int t) const;
	void dumpR_regularized(unsigned int t) const;
	void dumpS_regularized(unsigned int t) const;

	void dump_all_Phi() const;
	void dump_all_FPhi() const;
	void dump_all_G() const;
	void dump_all_FG() const;
	void dump_all_R() const;
	void dump_all_S() const;

	void dump_all_FPhi_regularized() const;
	void dump_all_R_regularized() const;
	void dump_all_S_regularized() const;

#endif

};
// *****************************************************************************************
#endif
